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O ' Abstract 



We consider short pulse propagation in nonlinear metamaterials characterized by a weak 
Kerr-type nonlinearity in their dielectric response. In the frequency "band gaps" (where 
linear electromagnetic waves are evanescent) with linear effective permittivity e < and 
^ I permeability > 0, we derive two short-pulse equations (SPEs) for the high- and low- 
Ph ' frequency band gaps. The structure of the solutions of the SPEs is also briefly discussed, and 
. connections with the soliton solutions of the nonlinear Schrodinger equation are presented. 
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, Artificially engineered metamaterials demonstrate unique electromagnetic (EM) prop- 

erties The frequency dependence of the effective permittivity e and permeability /i of 
lO ' these media is such that there exist frequency bands where the medium displays either a 
CSl ■ right-handed (RH) behavior (e > 0, /i > 0) or a left-handed (LH) behavior (e < 0, /i < 0), 
thus exhibiting negative refraction at microwave [I, 0, 0] or optical frequencies Fre- 
Q ! quency band gaps, i.e., frequency domains where linear EM waves are evanescent (e.g., for 
e < and > 0), also exist in metamaterials. Hence, when a nonlinearity occurs, say in the 
k>( ' dielectric response of the medium (a physically relevant situation in nonlinear metamateri- 
5_| ■ als [g], Q, [sI, M, E, 11, 12]), then nonlinearity- induced localization of EM waves is possible. 



Such localization is indicated by the formation of gap solitons, which occur mainly in non- 
linear optics and Bose-Einstein condensates (BECs) [l3], by means of the nonlinear 
Schrodinger (NLS) equation with aperiodic potential. Gap solitons were also predicted to 



occur in nonlinear metamaterials [15|. There, the approximation of slowly varying electric 
and magnetic field envelopes, led to a nonlinear Klein-Gordon (NKG) equation supporting 
gap solitons. 

Nonlinear models describing localization of wave packets in periodic media, e.g., the 



NLS equation in optics [13| and NKG equation in metamaterials [15| are usually derived in 
the framework of the slowly-varying envelope approximation. However, as far as ultra-short 
pulse propagation is concerned, i.e, for pulse widths of the order of a few cycles of the carrier 
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frequency, the NLS or the NKG models may faiL Indeed, in the context of nonhnear fiber 
optics, the proper model describing the evolution of such ultra-short pulses was shown to be 



the so-called short-pulse equation (SPE) [l6j. Numerical simulations in Ref. [16|| comparing 
solutions of Maxwell's equations to the ones of the SPE and NLS models have shown that the 
accuracy of the SPE (NLS) increases (decreases) as the pulse width shortens. More recently, 
motivated by the fact that the SPE model has no smooth pulse solutions propagating with 



fixed shape and speed, a regularized SPE (RSPE) was derived [l7|. The RSPE supports, 
under rather strict conditions, smooth traveling wave solutions. 

In this work, we derive nonlinear evolution equations describing ultra-short pulses that 
can be formed in the band gaps of nonlinear metamaterials. In doing so, we consider a 
metamaterial characterized by the permittivity e and permeability n of Ref. Ml, as well as 



a weak Kerr- type nonlinearity in the medium's dielectric response [10|, [Uj, |15|]. Following 
the lines of Refs. 16, we derive appropriate expressions for /i in the high-frequency 
(HP) and low-frequency (LP) band gaps with e < and /z > 0. Then, we use a multiscale 
perturbation method, with different small parameters for the HP and LP gaps depending 
on the metamaterial characteristics, to derive from Maxwell's equations two SPEs. Each 
of these equations describes the evolution of ultra-short pulses either in the HP or the LP 
gap and can not be valid simultaneously for a given metamaterial. We also discuss the 
structure of the resulting peakon-like solutions of these equations, and draw parallels to 
NLS-like soliton solutions (which can be regarded as ultra-short gap solitons in nonlinear 
metamaterials). Such peakon-like and breather-like solutions are also briefly considered 
numerically both in the context of the SPE models and in that of the originating Maxwell's 
equation dynamics. 

We consider a composite metamaterial consisting of an array of conducting wires and 
split-ring resonators (SRRs), characterized by a frequency dependent effective permittivity 
e{u!) and magnetic permeability fi^uj) (below we will use / and / to denote any function 
/ in the time- and frequency-domain, respectively). Here, the functions i^uj) and fi{uj) are 
assumed to take the form (in the lossless case) [2], 



e{uj) = eo 1 r , ij,{uj) = /io 1 
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where eo and fio are the vacuum permittivity and permeability, u is the frequency of the EM 
wave, ujp the plasma frequency, F the filling factor, and u^es the resonant SRR frequency. 
Assuming that ujp > Uj-^s, the dispersion relation = u^e{u)fi{u) {k is the wavenumber) 
shows that linear waves exist (fc^ > 0) in the bands oj > ojp and uj^^s < < ujm = 
uj^es/ Vl — F, where the medium displays right-handed (RH) and left-handed (LH) behavior, 
respectively. On the other hand, there exist two frequency bands, with e < and fi > 0, 
where linear EM waves are evanescent (fc^ < 0). These frequency "band-gaps" are < u < 
Clares ["low-frequency" (LP) gap], and um < uj < Up ["high-frequency" (HP) gap]. Typical 
dispersion curves, showing e{uj) and fi{uj), are depicted in Pig. [H 

We now consider a nonlinear metamaterial, which can be realized by filling the slits of 
the SRRs with a weakly nonlinear dielectric ll2|. In particular, we assume that 
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Figure 1: (Color online) Linear parts of the relative magnetic permeability, fi/iio [solid (red) line], and 
electric permittivity, e/eo [dashed (blue) line] as functions of the frequency uj in arbitrary units (a.u.). 
Propagation of linear waves is allowed (e/t > 0) for uj > tOp and Wres < to < u>m, where the medium displays 
right-handed (RH) and left-handed (LH) behavior, respectively. Dotted regions with e < and jj, > 0, 
namely, < lu < lo^cs (LF gap) and ujm < uj < cUp (HF gap), are frequency band gaps where linear waves 
are evanescent. 



this metamaterial exhibits a weak cubic (Kerr-type) nonhnearity in its dielectric response 



id, lUl, uM ' described by a nonhnear polarization vector Ptvl of the form, 

/ + 00 
XNL{t - ri,t - T2,t- r3)(E(ri) ■ E(r2))E(r3)rfrirfr2c/r3, (2) 

where E is the electric field, and xnl is the nonlinear electric susceptibility of the medium. 
In the case of small-amplitude, ultra-short pulse propagation, the nonlinear response is 
instantaneous, namely, 

XNL{t -ri,t- T2,t- ra) = K5{t - Ti)5{t - T2)5{t - Tg), (3) 

where n is the Kerr coefficient given hy K = ±^-^ with being a characteristic electric 
field value (e.g., of the order of 200 V/cm for n-InSb j§|); generally, both cases of focusing 
(k > 0) and defocusing (k < 0) dielectrics are possible. Notice that Eqs. ([2j) and ([3]) imply 
that Ptvl = eoK(E ■ E)E. 

We now assume propagation along the +z direction of a x- (y-) polarized electric (mag- 
netic) field, namely, E{z,t) = ^E{z,t) and li{z,t) = yH{z,t), where x and y are the 
unit vectors in the x and y directions, respectively. Then, Maxwell's equations lead to the 
following nonlinear wave equation for E{z,t): 

dlE -d'f{e*i2*E)- eoKd'^ijJ * E^) = 0, (4) 

where * denotes the convolution integral, f{t) * g{t) = f{T)g{t — T)dT, of any functions 
f{t) and g{t). Note that once the electric field E{z, t) is obtained from Eq. (jlj), the magnetic 
field H{z, t) can be derived from Faraday's law. 
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We first consider the high-frequency (HF) band gap, ujm < uj < Up, and assume EM 
wave frequencies u such that u ^ u^cs- In this regime, fi{uj) in Eq. ([1]) is approximated by: 

/iM^/io(l-F)-/ioi^4'. (5) 

Using the physically relevant parameter values [ol F = 0.4, Uj-es = 27r x 1.45 GHz, cum = 
2ti X 1.87 GHz, and = 271 x 10 GHz , in Fig. [2] we show the exact [Eq. ([T])] and 
approximate [Eq. ([5])] expressions for the effective permeability in this band. As seen, the 
above approximation produces a relative error from the exact form of (x{uj)/ixq less than 
5% in a wide sub-interval of frequencies in this band, namely for uj^, = 27r x 3.1 GHz 
< a; < cUp = 27r X 10 GHz. 
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Figure 2: (Color online) Linear part of the relative permeability /i//io in the HF gap {F = 0.4, Wies = 
27r X 1.45 GHz, ujm = 27r x 1.87 GHz, and = 27r x 10 GHz). Sohd (red) and dashed (blue) hnes 
correspond, respectively, to the exact [Eq. ([1])] and approximate [Eq. ©] expressions of in this 

band. The approximation produces a relative error less than 5% for Wa = 27r x 3.1 GHz < < o^p = 27r x 10 
GHz. 

The expression of Eq. ([5]) is useful for simplifying terms of Eq. (j4]) involving convolution 
integrals. In particular, the terms e* fi* E and fi* can respectively be approximated in 
the frequency domain as (l/c2)(l - F)E - {l/uj^c^)[{l - F)ujI + Fu'^^jE and {1/c^)k{1 - 

F)E^ — {1 / u'^c^) kFu'^^J^^ . Here, E = E exp^iujt) dt is the Fourier transform of E and 
c the velocity of light in vacuum. As a result, Eq. (jlj) can be reduced to the form: 

dlE - ^-^d^E - i [F< + (1 - F)ujI]E - ^ [Fu;LE' + (1 - F)d^E'] = 0. (6) 

Next, measuring time, space, and the field intensity E"^ in units of uj~l, v/oj^cs [where v = 
c(l — F)~^/^], and respectively, Eq. ([6]) is expressed in dimensionless form as 

{d!-d^-a)E = s(^-^ + d^E^ (7) 

where s = sgn(K) = ±1 for focusing or defocusing nonlinearity, respectively, and a = 

F/{1 — F) + {uopl LOres)"^ ■ Wc assumc that {ujp/uo^c&Y ^ 1 (as, e.g., in Ref. in order to 
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ensure the validity of the approximation ([5]) in a wide sub-interval of the HF band gap (see 
above). Hence, a is a large parameter, which suggests that a = a/e, where e is a formal 
small parameter (which sets also the field amplitude as per our perturbative approach below), 
and a = 0{1). Furthermore, considering propagation of small-amplitude short pulses, we 
introduce a multiple scale ansatz of the form 

E = e^/^E,{TuF, Z,,---)+ e'/^E2{TuF, Z,,- ■■) + ■■■ , (8) 

where Thf = s^^{t — z) and = e^z [n = 1,2, ■ ■ ■). Substituting Eq. (jS]) into Eq. ([7j), 
we obtain various equations at different orders of e. In particular, terms at 0(6'^^"^) cancel, 
there are no terms at 0{e~^^'^), while terms at 0{e~^^'^), cancel provided that the field Ei 
satisfies the following equation, 

2dcdT^^Ei + aEi + sdl^^El = 0, (9) 

where we have used the notation Q = Zi. Equation ([9]) is the so-called SPE, which was 



derived in Ref. [16| as an appropriate model describing the propagation of ultra-short pulses 
in silica optical fibers with a Kerr nonlinearity. 

Next, consider the low-frequency (LF) band gap, < to" < a;res, and assume that the EM 
frequency is a; <^ uj^es- Then, fi{uj) in Eq. ([1]) is approximated by 

/i(o.)^/xofl + F^y (10) 



Using F = 0.02, cjres = 27r x 1.45 GHz, in Fig. Owe show the exact [Eq. ([T])] and approximate 
[Eq. (II op ] expressions for fi{uj) in the LF band gap. This produces a relative error less than 
5% in a wide sub- interval of frequencies in this band, i.e., for < < cjb = 27r x 1.28 GHz. 
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Figure 3: (Color online) Linear part of the relative permeability fi/^io in the LF gap {F = 0.02 and 
Wres = 27r X 1.45 GHz). Solid (red) and dashed (blue) lines correspond, respectively, to the exact [Eq. ([1])], 
and approximate [Eq. (|10p ] expressions of jj,{u!)/fio hi this band. The approximation produces a relative 
error less than 5% for cjb = Stt x L28 GHz < a; < lu^cs = 2n x L45 GHz. 

We now employ Eq. flTUl) to simplify the convolution integral terms of Eq. (jl]). The 

terms e * fi * E and fi * E^ can respectively be approximated in the frequency domain as 
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(l/c2)[Fc^V^res + (1 " ^^p/^res) " ^p/^']^ and {1/c^)k{1 + Ft^V^^^jES. As a result, Eq. 
(HD reads, 

8lE - i (l - f ^) a?B - f E = + - ^,8tE\ (11) 

Notice that in Eq. (fTTi) . the ratio (cUp/c^res)^ is considered to be a C(l) parameter (as, e.g., 
in Ref. [3]) since ujp is not involved in the band width of the LP band gap. In this band, 
it is convenient to use a different sm all p arameter, namely the filling factor F, which is 



a physically relevant choice for SRRs [18|, [lOj, as well as for other types of metamaterials 



20, l2l|, 122] . Then, measuring time, space, and the field intensity E'^ in units of uj^^, c/uj^ 



and \k\ ^, respectively, we reduce Eq. f llll) to the dimensionless form: 



dl-dl-!^]E = sc)tE\ (12) 



v_ 

,2 

'res 

Next, we can again derive from Eq. 0121) a SPE for the LP band using the asymptotic 
expansion 

E = £Fi(Tlf, ^1, ■ ■ ■ ) + ^'^2(Tlf, Zi, ■■■) + ■■■ , (13) 

where Tlf = e^^{t — z) and Z„ = e'^z (n = 1, 2, ■ ■ ■ ). Then, substituting Eq. (|T3|) into Eq. 
(fT2|) . we find that terms at 0{e^^) cancel, there are no terms at C(£°), while terms at C(e), 
cancel provided that Ei satisfies the following SPE, 

2d^dT,,E, + ^E, + sdl^^El = 0, (14) 

^rcs 

where again ( = Zi. In the above analysis, the filling factor F was treated as 0{e^), with 
j > 5. However, if the filling factor F was assumed to be 0{e'^) then the additional term 
-d^^^Ei would appear in the left-hand side of Eq. (IT^ . In such case, Eq. (HM would then 



be the so-called regularized SPE (RSPE) model, which was recently derived in Ref. 17 



also in the context of ultra-short pulse propagation in nonlinear optical fibers. However 



due to its negative sign, the term —d^^^Ei does not have the regularizing effect of 17l |. 
and higher-order regularizations, outside the scope of the present work, may need to be 
considered. 

Solutions of the SPEs are now briefly discussed. Pirst, we unify SPEs of Eqs. ([9]) and 
(fT^ in the single equation 

d^dru + 7M + ^sd^^u^ = 0, (15) 

where u = Ei, while 7 = a/2 and r = The, or 7 = u!p/{2u!'^^^) and r = Tlf for the HP or 
the LP band gap, respectively. Then, seeking traveling wave solutions of the form u = u{^), 
where ^ = C ~ Ct (with C being associated with the velocity of the traveling wave), Eq. 
(ITB]) is reduced to the following ordinary differential equation: 



— Cu^^ + 'yu -\ — sC^(6'U'u| + 3u^u^^) = 0. 
^ 6 



(16) 



The transformation 
be solved to give 



w{u) produces a linear equation with respect to w{u) which can 



7U 



2 3sCu^ 



C (3sCm2-2) 



(17) 



subject to the initial condition ^^(±00) = m(±oo) = or ^(^(±00)) = w{0) = 0. The sign 
of the product sC is rather important. Indeed, if sC < then the solutions are always either 
ascending or descending since no maxima or minima can occur. Thus, we are left with the 
choice sC > 0, which allows bounded solutions. In the case s = +1 (i.e., focusing dielectrics 
with /€ > 0), which implies that C > 0, the maximum of the traveling wave occurs when 
= or, equivalently, u = ^J^j^C . To avoid the singularity at m = ^J2/?>C we consider 
small amplitude pulses (essentially pulses that never reach the singularity) and Eq. (ITTj) is 
reduced to the equation: 



{■j/C)u'^. The latter, possesses a peakon-like solution (see. 



e.g. 



Ref. 



of the form, 



«(0 = ^exp(-V77^iei), (18) 

whose derivative has a discontinuity at ^ = 0, with amplitude A < a/2/3C. Also, when the 
above approximation is not used, the right-hand side of Eq. 
a range of values for the field 



T7j) needs to be positive giving 
These observations are consistent with the 
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\u\ < y^I/Wl 

fact that the SPE exhibits loop-solitons found in Ref. 

To further illustrate, we now conduct a phase plane analysis of Eq. (|T6|) . see Fig. S] 
(C = s = 1). From the relevant curves, it becomes clear that there is no homoclinic orbit 
surrounding the fixed point at the origin. The only possibility is that we move on one of 
the phase plane curves, say, in the upper half plane up to a certain point, then "jump" 
from {u, U() to {u, —u^) and then, due to reversibility, return along the symmetric curve in 
the lower half plane. This would constitute the peakon-like solutions discussed above in the 
regime where u is small. However, we should also note that when integrating Eq. f|T5l) . we 
were not able to observe robust propagation of such a waveform in the dynamics (hence it 
is not discussed further herein). 




Figure 4: Typical phase plane curves {u,u^) associated with Eg. p^ for C = s = 1. 
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In addition, based on the formal connection between the SPE and the sine-Gordon 
equation (SGE), a smooth approximate, sech-shaped, envelope soliton solution of the SPE, 



based on the breather solution of the SGE, was derived in [2J]. In the framework of Eq. 
(fT5|) . this solution has the approximate form 



u 



4m(37)"^/^ cos (C + 7r) sech [m (C - 7r)] , (19) 



where m is an arbitrary real parameter, < m < 1. It is clear that the shape of the SPE 
pulse in Eq. ( IT9|) bears resemblance to the NLS soliton, as it consists of a sech-shaped 
pulse modulating a periodic function. The existence of this approximate envelope soliton 
solution, characterized by a small amplitude and inverse width (both determined by the 
parameter m < 1), suggests a connection between the SPE and the NLS equation. Such a 
connection can be established using the method of multiple scales as follows: introducing 
the variables Cn = ^"C) '^n = (with e <^ 1 and n = 0, 1, 2, ■ ■ ■ ) and expanding the field 
as M = '^n=i ^^"^niTn, Cn) , wc obtain from Eq. f|T5l) (for s = +1) the following results. The 
unknown field ui is found to be of the form 

Ml = A{(i, ri, ...) exp[i{k(o - lutq)] + c.c. (20) 

where c.c. denotes complex conjugate, while the wavenumber k and frequency u are con- 
nected by the dispersion relation [found to order C(e)]: uk + 'j = 0. On the other hand [as 
found at 0{e^)], the unknown envelope function A satisfies the following NLS equation, 

i {dc,A + k'dr.A) - ^d\A + ^uj\A\''A = 0, (21) 

where fi = Ti — k'd, k' = dk/du = —k/u, and k" = d'^k/du'^ = 2k/uj'^. Thus, it is clear 

that the well-known sech-shaped envelope soliton solution of the NLS Eq. (1^ resembles 

the soliton of Eq. f|T9l) . and scales in space and time in a similar way (if m is of 0(e)). Such 

smooth solutions of the SPE models derived above can be regarded as weak gap solitons (in 

the respective scales) that can be formed in the HP and LP band gaps of the considered 

nonlinear metamaterials. 

To corroborate these results, we have performed numerical simulations of both Eq. f[T^ . 

as well as of Eq. ([7]) from which the former was derived. Our numerical method relies 

on Fourier transforming Eq. (|T5l) with respect to r, then solving the ensuing first order 

ODE in ( (for each frequency), via a fourth-order Runge-Kutta scheme, and then Fourier 

transforming back to obtain u{(,t). In the simulations below, we use s = 1, F = 0.4, and 

^p/^res = 10/1.45 and e = 0.1. The initial condition is shown in the top left panel of Fig. \5\ 

and is obtained as the exact breather solution of the SPE equation (Eq. (22) in [24], with 

m = 0.32). The evolution of the breather can be seen both from the bottom left panel of Fig. 

[5l showing the center of mass of the solution vs. time and the top right panel illustrating 

the contour plot of the full space-time evolution. This is clearly a robust localized structure 

which propagates through the domain with constant speed in time (under the used periodic 

boundary conditions). We have also integrated Eq. ([7]) with the same type of breather-like 

initial profile. In the latter case, however, from the multiple scale ansatz, the initial condition 
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of Eq. ([7j) was chosen as E{0,t) = eu{0,T), EyjO, t) = —Ur{0,T). Furthermore, here one 
needs to be careful, similarly to what was done in |25| , to eliminate the propagation at very 
low frequencies [below \/& in the setting of Eq. ([7j)]. The result of the time integration 
is shown in the bottom right of Fig. [5] -see also the solid line showing the center of mass 
evolution in the bottom left panel of the figure. It is clear that the breather is robust in 
this setting as well, although its propagation speed is slightly smaller than that of the SPE 
breather. This result confirms our prediction that such "gap breathers" should be observable 
in nonlinear metamaterials of the type considered in this work. 
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Figure 5: The top left panel shows the breather initial condition used in Eqs. ^TE\i and The bottom left 
panel shows the evolution of the center of mass of the breathers in the two respective models by dotted and 
solid lines. The top right panel shows the space-time contour plot of the field evolution with the breather 
initial condition in Eq. (|15p . The bottom right panel shows the same for Eq. ([7]). 

In conclusion, we derived short-pulse equations (SPEs) describing the propagation of 
ultra-short pulses in nonlinear (Kerr-type) metamaterials. Two SPEs were found for the 
high- and low-frequency band gaps, respectively, characterized by a negative (positive) hn- 
ear effective permittivity (permeability), where propagation of linear electromagnetic (EM) 
waves is not allowed. We also discussed the structure of the solutions of the SPEs and 
presented the approximate peakon-like and breather-like solitary waves, which can be re- 
garded as weak ultra-short gap solitons. We also examined these structures via numerical 
computations to illustrate the apparent non-robustness of the former, and stable propaga- 
tion of the latter. Generally, the existence of such structures, indicates the possibility of 
nonlinear localization of EM waves in the gaps of nonlinear metamaterials. Interesting sub- 
jects for future research would include systematic studies of the stability of such ultra-short 
gap solitons both in the framework of the SPEs and Maxwell's equations as well as higher 
dimensional generalizations of the structures considered herein. 
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